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ABSTRACT 


Conputar  aidad  daaign  of  nonllnaar  dlfferantlal  ayacens  require 
■any  raaponaa  calculatlona  from  the  ayatam  modal  x  ■  F(x,u).  Many 
nonllnaar  ayatama,  aapaclally  control  ayatama,  can  ba  aaparated  Into 
a  atiff  llnaar  part  and  a  nonatlff  nonlinear  part.  Faat  numerical 
Integration  tachnlquea  taking  advantage  of  the  partitioned  form  are 
praaantad. 


I.  INTRODUCTION 


Computer  aided  design  of  nonlinear  differential  systems  requires 
many  response  calculations  from  the  system  model  x  ■  F(x,u)  where 
X  Is  the  state  vector  and  u  Is  the  control  vector.  Accurate  models 
Include  both  short  and  long  term  effects,  and  thus  have  widely  separated 
or  stiff  eigenvalues.  For  large  stiff  sets  of  equations  the  analysis 
takes  excessive  execution  time  as  most  Integration  techniques 
(Including  the  Runge  Kutta  rules)  require  that  the  integration  time 
step  be  limited  by  the  magnitude  of  the  largest  and  possibly  most 
uninteresting  eigenvalues. 

Many  nonlinear  differential  systems,  especially  control  systems, 
can  be  separated  Into  a  linear  and  a  nonlinear  part  of  the  form 
X  ■  Ax  +  f(x,u)  where  all  the  stiff  eigenvalues  are  In  the  linear 
part  represented  by  the  matrix  A.  This  partitioning  is  possible 
for  most  cases  because  the  nonlinear  equations  usually  describe  the 
motion  of  the  system  to  be  controlled  (e.g.  airplane,  missile,  hydro¬ 
foil,  etc.)  which  has  time  responses  much  slower  than  those  of  the 
controlling  sensors,  compensators  and  actuators.  Further,  most 
sensors,  compensators  and  fast  actuators  are  designed  to  operate  In 
a  linear  or  almost  linear  region  of  their  response  curves.  Non- 
llnearltles,  such  as  physical  limits  on  actuator  travel,  can  most 
often  be  Incorporated  Into  the  plant  equation  as  Input  saturation. 

Thus  the  total  system  Including  feedback  which  Is  represented  by 
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F(x,u)  can  be  partitioned  so  that  the  large  eigenvalues  occur  in  the 
linear  part. 

In  the  numerical  response  calculation  for  nonlinear  differential 
systems,  there  are  three  types  of  errors  to  consider;  truncation  error, 
roundoff  error,  and  error  due  to  numerical  instability  of  the  integrating 
rule.  Roundoff  error  can  be  controlled  by  use  of  sufficient  word  length 
and  programning  techniques  such  as  inner-product  accumulation.  Truncation 
error  is  the  one  step  error  due  to  the  diacrete  approximation  to  the 
Integral.  This  error  is  a  function  of  the  integration  rule,  the  time 
step  size,  and  the  state  of  the  ayatem.  For  a  particular  integration 
rule  and  step  size,  bounds  for  the  truncation  error  can  be  determined, 
however,  actual  truncation  error  can  be  significantly  less  depending 
on  the  state  of  the  system.  For  example,  if  the  state  of  the  system 
is  a  stable  singular  point  of  the  nonlinear  equations,  truncation 
error  could  be  negligible.  The  third  type  of  error,  that  arising 
from  numerical  Instability  of  the  integration  rule,  is  often  dis¬ 
missed  because  truncation  error  control  schemes  will  normally  force 
the  step  size  small  enough  so  that  the  integration  rule  is  stable. 

For  stiff  systems  this  nay  cause  an  unreasonably  large  nuri>er  of  inte¬ 
gration  steps  to  be  required  by  some  rules. 

Numerical  stability  of  integration  rules  [1,2]  normally  depends 
on  the  product  of  the  maximum  modulus  eigenvalue  (of  the  Jacobian)  and 
the  time  step  size  T.  For  exaaq)le,  for  a  Runge  Kutta  fourth  order  rule 
|x  max]  T  must  be  less  than  2.8.  Thus  even  though  the  truncation  error 
might  be  negligible  for  larger  step  sizes,  T  must  remain  ssull  to 
guarantee  the  computed  solution  does  not  grow  without  bound. 

The  purpose  cf  this  paper  is  to  present  several  integration  rules 


whlch  are  stable  over  much  larger  regions  than  standard  rules  and 
thus  allow  the  step  size  to  be  Increased  when  the  truncation  error 
allows.  For  these  rules,  the  region  of  numerical  stability  depends 
on  the  maxiaum  modulus  eigenvalue  of  the  nonlinear  part  only,  which 
due  to  partitioning  may  be  orders  of  magnitude  less  than  that  of  the 
composit  system. 


II.  RESPONSE  CALCULATION 


Transformation 

Assum  a  system  of  the  form  x  -  Ax  +  f(x,u)  where  stiff 
eiganvaluwS  are  only  in  the  linear  part  A.  Transforming  the  state 
variables,  x  ■  e^S*  gives  a  new  set  of  nonstiff  differential 
equations,  y  •  G(y,u)  obtained  as  follows, 


X  ■  Ax  +  f(x,u) 

,  At  .  At*  *  At  .  At  V 
Ae  y  +  e  y«Ae  y  +  f(e  y,u) 

•  -At-.  At  V 

y  ■  e  f(e  y,u)  . 


The  dynamic  eigenvalues  of  a  nonlinear  system  are  the  eigen¬ 
values  of  Che  Instantaneous  Jacobian.  The  eigenvalues  of  the 
original  system  are  the  eigenvalues  of  A  f^(x,u).  The  eigen¬ 
values  of  the  transformed  system  are  the  eigenvalues  of 

Gy  -  e  ■  ®  ^x’^y  *  f^{x,u)e  (4) 

Since  G  Is  a  similarity  transformation  of  f  ,  the  eigenvalues  of 

y  X 

the  transformed  system,  G,  are  the  eigenvalues  of  the  nonlinear  part 


I 


of  the  original  system,  £. 


The  transformed  system  of  equations  can  be  integrated  by  fast 

explicit  integration  techniques  without  requiring  small  time  steps 

to  prevent  numerical  instabilities  due  to  stiff  eigenvalues.  Inte- 

grating  y  over  one  time  step,  t  . ,  -  t  ■  T,  gives 

nTi  n 


/nfl 


Premultiplying  by  e 


A^n+l 


gives  the  convolution  integral  sequential 


difference  approximation  to  x^^^ 


AT 

X  .  1  ■  e  X 
n+1  1 


f(x,u)dt 


Numerical  stability  is  maintained  by  using  either  an  exact  analytic 
AT 

expression  for  e  or  by  using  the  numerically  stable  Fade  approxi- 

AT 

matlons.  Table  1  presents  some  symmetric  Fade  approximations  to  e  . 

AT 

Numerical  techniques  for  the  analytical  determination  of  e  using 
similarity  transformations  are  similar  to  those  presented  in  [3]. 


Table  1.  Numerically  stable  symmetric  Fade 


approximations  to  e  . 


Order 


2  (I  -  yTA)"^  (I  +  "ItA) 

4  (I  -  -jTA  +  (I  +  ^A  + 

6  (I  -  |tA  +  ^tV  -  (I  +  |tA  +  4.  ^tV) 

8  (I  -  ^TA  +  IgTV  -  I^tV  +  +  |tA  4-  I^T^A^ 


4-  — T^A^  4-  — ^ — T^A^) 
84  1680  ^  '' 
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then  the  truncation  error  for  a  p-th  order  nile  using  a  p-  :h  order 
Fade  approximation  Is  dominated  by 

e  (11) 

n  n 

AT 

since  I |a| I  >>  ||b||.  If  the  approximation  to  e  Is  of  order 
greater  than  p  then 

e  -  KT^^^a’^Bx  .  (12) 

n  n 

For  T  much  greater  than  |  |Aj  |  ^  the  error  gets  large  so  that 
significant  Increases  In  T  cannot  he  made. 

For  example  if  A  haf<  eigenvalues  like  'lODO  and  f^  has 
eigenvalues  like  -1  then  a  standard  Runge  Kutta  rule  will  require  a 
time  step  of  about  .001  to  maintain  numerical  stability.  The  trans> 
formed  method  will  be  numerically  stable  for  any  step  slse  up  to  1» 
but  the  truncation  error  can  be  like  T***^10^**/(p+l) ! .  To  control 
this  truncation  error  T  might  still  have  to  be  kept  small. 

Convolution  Integral  Partitioning 

When  the  stiff  system  of  equations  can  be  further  partitioned 
so  that  part  of  the  system  has  no  stiff  part  as  Js  shown  In  1 

and  equation  13,  then  low  error  convolution  Integral  approximations 
can  be  made. 
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Flgure  1.  DIAGRAM  OF  PARTITIONED  SYSTEM 
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Nuacrlcally  acable,  low  error  intetretion  rules  can  be  obtained  by 
applying  first  and  second  oriUer  explicit  one  step  end  point  approx¬ 
imations  to  the  convolution  Integral  representation  of  X2> 

Vl 

x^l  -  e^2Tj^  J  f2(t)dt  (14) 

t 

n 

Tills  approxlaatlon  to  x^  Is  coupled  with  cosipatlble  Range  Kutta 
rules  to  approxlsMte  x^. 

A  first  order  approxlmstlon  to  the  convolution  Integral  can  be 
obtained  by  expanding  f^  about  t^  and  then  Integrating  the  ex¬ 
ponential  tern  exactly. 

f,(x(t))  -  fV  0(T)  (15) 

.  1  .  T 

J  .  I)  .  (16) 

t 

n 


..n>  ,■  ^ 
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The  integral  1^  is  well  defined  even  if  A2  is  not  of  full 
rank  and  can  be  calculated,  as  before,  using  stable  Fade  approxi¬ 
mations  or  evaluating  exact  analytical  expressions.  Using 
equations  (15)  and  (16),  equation  (14)  becomes 

+  O(T^)  .  (17) 

From  a  first  order  rule 

-  x"  +  TfJ  (18) 

Similarly  to  get  a  second  order  rule 

f2(x(t))  -  ^2  ^2  ^ 

y"  ^  -  TI)  -  I2  (20) 

t 

n 

A  first  order  approximation  to  f,  is  needed  so  evaluating  f^  at 

from  (17)  and  (18)  gives  a  first  order  approximation 

f^^  -  f 2  * 

T  “ 

(19),  (20),  and  (21)  equation  (14)  becomes 


n+1 


to 


and  X 
^n+l 


n>l 


Using  equations 


n+l 


K.om  a  second  order  Runge  Kutta  rule 


‘1  ■  *1  *  2^<'l  *  '1  > 


(23) 


|»  mim 
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where  again  is  evaluated  at  the  first  order  points  in 

equations  (17)  and  (18).  Equations  (22)  and  (23)  define  a 
partitioned  second  order  rule. 

Under  the  special  circumstances  when  is  a  function  of 

only  and  the  Jacobian  of  f^  is  readily  available  (such  as 

when  f^  is  linear)  then  slightly  better  accuracy  can  be  obtained  by 
using  the  following  approximation  to 

f,(x.(t))  -  f^  +  f"  f"  [t-t  ]  +  O(T^)  (24) 

Z  1  1  z\  L  n 


Now  equation  (22)  takes  the  form 


n+1  A9T  n  .  T  ,n  .  -  ,n  -n 
'2  -  »  2  Xj  +  Ijfj  +  1/2^  f, 


(25) 


Using  (18)  and  (25)  f^^^  can  be  approximated  to  be  used  to  define 

x^^^  in  equation  (23).  Equations  (22)  and  (23)  define  a  partitioned 
second  order  rule  when  the  Jacobian  of  is  available. 

The  truncation  errors  for  these  partitioned  rules  are  on  the 
order  of  the  truncation  errors  for  a  system  of  equations  with  eigen¬ 
values  on  the  order  of  the  eigenvalues  of  tlie  nonlinear  parts.  For 
example  if  A2  has  eigenvalues  like  -1000  and  f^  and  have 

eigenvalues  like  -1  then  tlie  truncation  error  will  be  like  (p+1)  ! 

To  get  3  place  accuracy  with  a  first  order  rule  T  should  be  about 
.04  and  with  a  second  order  rule  T  should  be  about  .1. 


III.  EXAMPLE 

In  order  to  illustrate  the  application  of  these  solution 
techniques,  consider  a  simple  pendulum  problem  in  which  the  angular 


> 


10 


rat*  i*  m*aaur«d  and  u**d  aa  a  algnal  to  a  torqua  notor  to  provld* 
danplng.  Ttia  nonllnaar  aquation*  daacrlblng  tha  notion  of  tha 
pendulun  ar* 

ML  0(t)  -  -Mg  ain  6(t)  ^  T(t)/L  (26) 


whare  M  ia  the  naaa,  L  ia  the  length,  and  T(t)  ia  the  torque 
applied  by  the  torqua  motor.  Let  the  dynamic*  of  the  angle  sensor 
be  given  by  a  first  order  lag  with  a  time  constant  Likewise 

let  the  torque  motor  be  represented  by  a  first  order  lag  with  time 
constant  1/t  .  To  obtain  the  state  equations,  let  x,  be  6, 

X2  ba  6,  be  the  output  of  the  angular  rate  sensor,  and  let 

be  the  torque.  The  equations  of  motion  are  then 


*1 

■*2 

X2 

-g/L  sin  Xj  -  x^/ML^ 

*3 

■ 

K  T  X.,  -  T„X, 
s  s  2  8  3 

K  T  X-  -  T  X. 

L  m  m  3  m  4  J 

(27) 


2 

For  the  purpose  of  this  example,  let  g/L  be  1.,  1/ML  »  0.5, 

T  ■  1000.,  K  -  1.0,  T  ■  50.,  and  K  ■  3.0. 

8  ’  8  ’  m  ra 

Now  factor  the  system  into  the  form 


X  ■ 


1 

- 

-X2 

- 

u 

0 

X  + 

sin  x^ 

-.5x4 

L  !  150 

- 1 

00 

m 

1 

CM 

X 

0  0 
0 

0 

_ 1 

(28) 


Ax  +  F. 
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Note  that  the  Jacobian  of 


F  for  small 

0  o' 

0  -.5 

0  0 

0  0 


(29) 


which  has  eigenvalues  of  0.0,  0.0,  ±J1.0.  The  eigenvalues  of  the 
composlt  system  for  small  are  -1000.08,  -48.371,  -.77484±j  .65818. 

For  the  example  system  of  equations  considered  here  some  Indi¬ 
cations  of  the  acceptability  of  the  suggested  Integration  rules  can  be 
made.  Figure  2  Is  a  plot  of  discretization  error  curves  for  several 
rules.  The  log  of  the  error  versus  the  log  of  the  Integration  time 
step  has  a  slope  which  Is  equal  to  the  order  of  the  Integration  rule. 
Table  2  compares  the  number  of  significant  figures  of  accuracy  obtained 
for  the  Integrating  rules  at  several  time  steps.  Also  computer 
execution  times  are  given. 

A  fourth  order  Runge  Kutta  rule  had  high  accuracy  for  small  step 
sizes,  but  It  went  unstable  for  acceptable  step  sizes.  An  order  two 
Runge  Kutta  rule  went  unstable  for  smaller  step  sizes.  The  trans¬ 
formed  Runge  Kutta  rules  remained  stable  for  large  step  sizes  but  the 
accuracy  was  poor.  The  partitioned  methods  ^ave  acceptable  accuracy 
for  reasonable  step  sizes.  Using  the  Jacobian  technique  gave  about 

1  place  better  accuracy  than  the  time  derivative  approximation  technique. 

AT 

Using  a  Fade  approximation  to  e  instead  of  an  exact  expression  did 
not  degrade  the  accuracy  much,  but  In  some  problestf  this  may  cause  loss 
of  accuracy  because  large  negative  Fade  exponents  approach  -1.  Instead 
of  0.  The  partitioned  2nd  order  rule  allows  integration  to  acceptable 
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accuracy  about  200  times  faster  than  a  standard  4th  order  Runge 
Kutta  rule  for  this  example. 


Table  2.  Number  of  significant  figures 
and  timings  for  example. 


Integration  rule/time  step 

1 

O 

10"^ 

Order  4  Runge  Kutta 

12 

unstable 

Order  2  Runge  Kutta 

unstable 

unstable 

Order  2  transformed 

3 

1 

Order  2  partitioned 

6 

4 

Order  2  partitioned 
(with  Fade  approximation) 

6 

4 

Order  2  partitioned 
(with  Jacobian) 

7 

5 

Execution  time  in  seconds 
for  all  order  2  rules 
(360/44) 

30 

3 

unstable 

unstable 


0 

2 

2 


3 


.3 
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Figure  2.  INTEGRATION  ERROR  CURVES 
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IV.  CONCLUSIONS 

In  many  simulations  vhere  one  or  .1  percent  accuracy  is 
sufficient,  the  larger  step  sizes  allowed  by  the  low  order  stable 
integration  rules  reduces  run  times  appreciably.  Control  systems 
for  nonlinear  plants  which  use  position  and  rate  feedback  with 
compensation  usually  allow  partitioning  as  indicated  in  Eq.  13  with 
f^  being  a  function  of  and  u  only.  Experience  has  shown  the 

second  order  rule  to  be  ideally  suited  for  this  problem.  Vfhere 
acceleration  feedback  is  used,  f^  now  becoMs  a  function  of  X2 
in  addition  to  x^  and  u  which  in  some  cases  reduces  the  effective¬ 
ness  of  the  integration  rule. 

The  rules  presented  here  show  advantages  over  other  rules  such 
as  standard  Runge  Kutta  only  when  the  truncation  errors  of  the  p*‘imary 
variables  are  small  enough  to  warrant  larger  step  sizes.  Vfhere  the 
problem  is  truncation  error  limited,  the  higher  order  rules  offer 
advantages  in  error  control.  It  should  be  noted  that  in  many  simulations, 
however,  errors  due  to  numerical  inst.sbillties  are  the  limiting  factor. 

For  these  systems,  stable  integration  rules  can  significantly  reduce 
computation  times. 

One  obvious  example  where  these  rules  are  particularly  effective 
is  in  determining  steady  state  operating  conditions  of  nonlinear  systems. 
Here,  large  errors  in  the  transient  solution  can  be  tolerated  sulking 
stable  integration  rules  with  large  step  sizes  quite  desirable. 


I 
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